function h = adf_lag(y , hmax , trend , plevel)
% Select number of lags h to use in the augmented Dickey-Fuller test:
%  y_t = \mu + \beta t + \pi y_{t-1} + \sum_{i=1}^h\gamma_i\Delta y_{t-i} + u_t
% Following Campbell and Perron (1991) and Hall (1994) Section 1.2
% Currently for univariate time series only
% INPUT:
%   y -- univariate time series being tested
%   hmax -- maximum number of lags to be considered
%   trend -- trend to be used in detrending. Default constant
%            1: include both a constant and a time trend
%            0: include only a constant
%            -1: no constant or time trend
%   plevel -- significance level. Default 0.05
%% Input Checking
if nargin < 2
    error('2 inputs required')
end
% kmax needs to be at least 1
if hmax < 1
    error('kmax needs to be at least 1')
end
% Get the length of the data
[T,K]=size(y);
if K > 1
   error('DATA must be univariate.')
end 
if T < hmax + 1 
    error('DATA does not have enough observations.')
end
% Default option for trend
switch nargin
    case 2
        trend = 0;
        plevel = 0.05;
    case 3
        plevel = 0.05;
end

%% Detrend
% RHS constant and trend term
if trend == 1
    RHS_exo = [ones(T , 1) , (1 : T)'];
    b_trend = olsgmm(y , RHS_exo , 0 , -1);
    ty = y - RHS_exo * b_trend;% This is the \tilde{y} in Campell and Perron (1991)
elseif trend == 0
    RHS_exo = ones(T , 1);
    b_trend = olsgmm(y , RHS_exo , 0 , -1);
    ty = y - RHS_exo * b_trend;
elseif trend == -1
    ty = y;
end

%% Estimate AR significance
% Get AR difference
ty_dif = ty(2 : end) - ty(1 : end - 1);
% Get lagged matrix for RHS
ty_dif_lag = nan(T - 1, hmax);
for j = 1 : hmax
    ty_dif_lag(j + 1 : end , j) = ty_dif(1 : end - j);
end
% Testing significance for each lag
reg_result = nan(hmax , 1);
for i = 1 : hmax
    mdl_temp = fitlm([ty(i + 1 : end - 1) , ty_dif_lag(i + 1 : end , 1 : i)], ty_dif(i + 1 : end) , 'Intercept',false);
    reg_result(i , 1) = mdl_temp.Coefficients.pValue(end);%store pvalue
end
reg_sig = reg_result < plevel;
%% Determine lags
% We use the specific-to-general approach in Campell and Perron (1991)
% If last included lag significant, k = kmax
% Otherwise include the last significant lag
h = find(reg_sig , 1 ,'last');
if isempty(h)
    h = 0;
end

end
